Bowfin Migration and Home Range Analysis
Greg Jacobs
21 FEB 2024
t0 <- Sys.time()
set.seed(2) # make workflow repeatable
library(tidyverse)
library(ggplot2)
library(ggExtra)
library(ggpubr)
library(kableExtra)
library(tidyr)
library(readxl)
library(lubridate)
library(sf)
library(raster)
library(amt)
library(purrr)
library(runjags)
library(patchwork)
library(lme4)
library(lmerTest)
library(emmeans)
library(sjPlot)
library(sjmisc)
library(sjlabelled)
library(DHARMa)
select<-dplyr::select
crs_utm18n <- "EPSG:26918" # UTM zone 18N, https://epsg.io/32618
lu <- function(x) length(unique(x)) # a helper function to save space
# fish metadata recorded after tagging and release
fishids <- read_excel("BowfinInventory_gj.xls") %>%
mutate(id=as.character(FishID))
# fish observation data
fishobs_tab <- read_csv("bowfin_data.csv") %>%
left_join(fishids) %>%
filter(!is.na(Latitude)) %>% # filter out missing coordinates
mutate(DOY=yday(DateTime)) %>% # need DOY of all datetime observations
filter(!is.na(distday)&Distance>0) %>% # keep only observations of movement
mutate(Sex=factor(Sex), Year=factor(Year), ReleaseSite=factor(ReleaseSite),
Zdoy=scale(.$endDOY)[,1])
fishobs <- fishobs_tab %>% filter(!is.na(Latitude)) %>%
st_as_sf(coords=c("Longitude", "Latitude"), crs=4326) %>%
st_transform(crs=crs_utm18n)
# Bayesian change-point models and results
load(file="mcp_out.RData")
# Bayesian change-point results summery table
mcp_summary <- read_excel("MCPresultstab_20240303.xlsx") %>%
select(-Notes) %>%
left_join(fishids)
# sf object of initial capture locations for each fish
caploc <- fishids %>%
mutate(x=if_else(ReleaseSite=="Bradburys", -76.11321, -75.92999),
y=if_else(ReleaseSite=="Bradburys", 43.25366, 43.17359)) %>%
st_as_sf(coords=c("x", "y"), crs=4326) %>%
st_transform(crs=crs_utm18n) %>%
mutate(Year=year(ReleaseDate), DOY=yday(ReleaseDate)) %>%
select(FishID, Year, DOY) Supporting spatial data
Load spatial datasets and project everything into NAD83/UTM zone 18N.
# Oneida Lake boundary from NHD v2.1
olb <- sf::st_read(dsn="supporting", layer="OneidaLake_utm18n") ## Reading layer `OneidaLake_utm18n' from data source
## `C:\Users\gj93\Box\Bowfin_telemetry\bowfin_analysis_mobile\bowfin_HR_SignerWorkflow_220709\manuscript_version\supporting'
## using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 14 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 407400 ymin: 4777457 xmax: 440730.1 ymax: 4789914
## Projected CRS: NAD83 / UTM zone 18Nlrwpoly <- sf::st_read(dsn="supporting",
layer="lrw_poly",
crs=crs_utm18n) ## Reading layer `lrw_poly' from data source
## `C:\Users\gj93\Box\Bowfin_telemetry\bowfin_analysis_mobile\bowfin_HR_SignerWorkflow_220709\manuscript_version\supporting'
## using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 2 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 400003.9 ymin: 4773813 xmax: 444803.9 ymax: 4793563
## Projected CRS: NAD83 / UTM zone 18Nbathy <- raster("supporting/depth_2")
rw <- raster("supporting/rwindicator")
lk <- raster("supporting/oneidaind")
lrw <- raster("supporting/lrwindicator")
# create data frames from rasters, clipped to lake shore, for plotting
bathy_df <- as.data.frame(bathy*lk, xy=TRUE)
names(bathy_df)[3] <- "value"
rw_df <- as.data.frame(rw*is.na(lk), xy=TRUE) %>%
mutate(value=ifelse(layer==1, 1, NA))
lrw_df <- as.data.frame(lrw*!is.na(lrw), xy=TRUE) %>%
mutate(value=ifelse(layer==1, 1, NA))
lrw_mask_df <- mutate(lrw_df, value=ifelse(is.na(value),1,NA))Map locations across available habitats in Oneida Lake and surrounding wetland and stream habitats.
# plot locations by FishID - base plot for fig 1 and results figures
plotlocations <- ggplot() +
geom_tile(data=rw_df %>% filter(!is.na(value)),
aes(x=x, y=y), fill='gray90', color='gray90') +
geom_tile(data=bathy_df %>% filter(!is.na(value)),
aes(x=x, y=y, fill=value, color=value)) +
scale_fill_viridis_c(name='Depth (m)', direction=-1, end=0.6) +
scale_color_viridis_c(name='Depth (m)', direction=-1, end=0.6) +
geom_sf(data=fishobs, shape=3, alpha=0.5) +
geom_sf(data=caploc, shape=19, col="gold") +
coord_sf(xlim=c(408000, 440000), ylim=c(4773813, 4793563)) +
theme_void() + #theme(legend.position="none")+
ggspatial::annotation_north_arrow(which_north="true", pad_y=unit(0.25, "in"),
style = ggspatial::north_arrow_fancy_orienteering) +
ggspatial::annotation_scale(location = "bl", width_hint = 0.5) +
theme(plot.title=element_text(hjust=0.5))Depth and Distance from shore.
Visualize and analyze sex differences.
fishobs_tab %>%
group_by(FishID, Sex) %>%
summarize(ShoreDist_m=mean(ShoreDist_m, na.rm=TRUE),
Depth_m=mean(Depth_m, na.rm=TRUE)) %>%
ungroup() %>%
group_by(Sex) %>%
summarize(mean_ShoreDist_m=mean(ShoreDist_m, na.rm=TRUE),
sd_ShoreDist_m=sd(ShoreDist_m, na.rm=TRUE),
n_ShoreDist_m=sum(!is.na(ShoreDist_m)),
mean_Depth_m=mean(Depth_m, na.rm=TRUE),
sd_Depth_m=sd(Depth_m, na.rm=TRUE),
n_Depth_m=sum(!is.na(Depth_m)))## `summarise()` has grouped output by 'FishID'. You can override using the
## `.groups` argument.## # A tibble: 2 × 7
## Sex mean_ShoreDist_m sd_ShoreDist_m n_ShoreDist_m mean_Depth_m sd_Depth_m
## <fct> <dbl> <dbl> <int> <dbl> <dbl>
## 1 F 130. 48.1 18 1.58 0.371
## 2 M 109. 63.0 18 1.35 0.227
## # ℹ 1 more variable: n_Depth_m <int>Repeated measures analysis of depth by sex.
lmer(log(Depth_m) ~ Sex + (1|FishID), data=fishobs_tab) %>%
summary(ddf="Satterthwaite")## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(Depth_m) ~ Sex + (1 | FishID)
## Data: fishobs_tab
##
## REML criterion at convergence: 1466.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.6265 -0.6147 -0.0067 0.6471 2.6218
##
## Random effects:
## Groups Name Variance Std.Dev.
## FishID (Intercept) 0.01494 0.1222
## Residual 0.30774 0.5547
## Number of obs: 864, groups: FishID, 36
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.28729 0.04372 34.68911 6.570 1.44e-07 ***
## SexM -0.14765 0.05857 31.17784 -2.521 0.017 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## SexM -0.747lmer(log(Depth_m) ~ Sex + (1|FishID), data=fishobs_tab) %>%
anova(ddf="Satterthwaite")## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Sex 1.9558 1.9558 1 31.178 6.3553 0.01703 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Repeated measures analysis of distance from shore by sex.
lmer(log(ShoreDist_m+1) ~ Sex + (1|FishID), data=fishobs_tab) %>%
summary(ddf="Satterthwaite")## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(ShoreDist_m + 1) ~ Sex + (1 | FishID)
## Data: fishobs_tab
##
## REML criterion at convergence: 2686.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.7910 -0.4359 0.1155 0.6736 2.2524
##
## Random effects:
## Groups Name Variance Std.Dev.
## FishID (Intercept) 0.1727 0.4155
## Residual 1.4405 1.2002
## Number of obs: 824, groups: FishID, 36
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 4.4551 0.1251 33.9960 35.619 < 2e-16 ***
## SexM -0.4843 0.1698 31.3363 -2.853 0.00761 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## SexM -0.737lmer(log(ShoreDist_m+1) ~ Sex + (1|FishID), data=fishobs_tab) %>%
anova(ddf="Satterthwaite")## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Sex 11.722 11.722 1 31.336 8.1375 0.007614 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1depthdistmvmntplot <- fishobs_tab %>%
group_by(FishID, Sex) %>%
summarize(barShoreDist=mean(ShoreDist_m+1, na.rm=TRUE),
barDepth_m=mean(Depth_m, na.rm=TRUE),
barDistDay=mean(distday, na.rm=TRUE))## `summarise()` has grouped output by 'FishID'. You can override using the
## `.groups` argument.pShoreDist <- depthdistmvmntplot %>%
ggplot(aes(barShoreDist, color=factor(Sex))) +
scale_color_grey() +
geom_density(size=1) +
scale_x_log10() +
ylab("Density") +
xlab("Shore Distance (m)") +
theme_classic() + theme(legend.position = "none")## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.pDepth <- depthdistmvmntplot %>%
ggplot(aes(barDepth_m, color=factor(Sex))) +
scale_color_grey() +
geom_density(size=1) +
scale_x_log10() +
ylab("") +
xlab("Depth (m)") +
theme_classic() + theme(legend.position = "none")
pdistday <- depthdistmvmntplot %>%
ggplot(aes(barDistDay, color=factor(Sex))) +
scale_color_grey() +
geom_density(size=1) +
scale_x_log10(breaks=c(1,10,100,1000)) +
ylab("") +
xlab("Movement (m/day)") +
theme_classic() + labs(color="Sex")
fig3 <- pShoreDist + pDepth + pdistday +
plot_annotation(tag_levels=list(c("a)", "b)", "c)"))) +
plot_layout(widths=c(1,1,1), ncol=3, nrow=1)
ggsave("fig3_freq.png", fig3, width=6.5, height=2, units="in", dpi=600,
bg="white")
knitr::include_graphics("fig3_freq.png")
Movement analysis
Analysis of minimum distance moved per day by sex, year, and day of year.
# population average and sd of individual average movement rate
mMdat <- filter(fishobs_tab, dayssince<22)
group_by(mMdat, FishID, Sex) %>%
summarize(meandistday=mean(distday, na.rm=TRUE)) %>%
group_by(Sex) %>%
summarize(meanM=mean(meandistday),
sdM=sd(meandistday))## `summarise()` has grouped output by 'FishID'. You can override using the
## `.groups` argument.## # A tibble: 2 × 3
## Sex meanM sdM
## <fct> <dbl> <dbl>
## 1 F 221. 157.
## 2 M 114. 107.# population average of within-individual SD of movement rate
mMdat <- filter(fishobs_tab, dayssince<22)
group_by(mMdat, FishID, Sex) %>%
summarize(sddistday=sd(distday, na.rm=TRUE)) %>%
group_by(Sex) %>%
summarize(meansdM=mean(sddistday, na.rm=TRUE))## `summarise()` has grouped output by 'FishID'. You can override using the
## `.groups` argument.## # A tibble: 2 × 2
## Sex meansdM
## <fct> <dbl>
## 1 F 280.
## 2 M 253.mM <- lmer(log(distday) ~ Sex*Year + Sex*Zdoy + (Zdoy|FishID), data=mMdat)anova(mM, ddf="Satterthwaite")## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Sex 64.961 64.961 1 34.91 21.3413 5.057e-05 ***
## Year 7.585 3.793 2 744.60 1.2459 0.288272
## Zdoy 8.989 8.989 1 31.25 2.9531 0.095606 .
## Sex:Year 30.519 15.260 2 744.60 5.0131 0.006876 **
## Sex:Zdoy 12.960 12.960 1 31.25 4.2576 0.047465 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1margMdf <- emmip(mM, ~ Sex|Year|Zdoy, lmer.df="satterthwaite", plotit=FALSE)
pd <- position_dodge(.3)
p <- ggplot(data=margMdf, aes(x=Year, y=exp(yvar), fill=Sex, Sroup=Sex)) +
geom_errorbar(aes(ymin=exp(yvar-SE), ymax=exp(yvar+SE)),
color='black', width=.2, position=pd) +
geom_line(position=pd) +
geom_point(color='black', size=8, position=pd, shape=21) +
scale_fill_grey() +
ylab(expression("Movement "~(m/day))) +
theme_classic()
ggsave("fig4_Mmeans.png", p, width=6.5, height=4, units="in", dpi=600,
bg="white")## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?knitr::include_graphics("fig4_Mmeans.png")
tab_model(mM, df.method = "satterthwaite")| log(distday) | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | 4.60 | 3.92 – 5.28 | <0.001 |
| Sex [M] | -2.09 | -2.98 – -1.20 | <0.001 |
| Year [2010] | -0.60 | -1.27 – 0.07 | 0.078 |
| Year [2011] | -1.02 | -1.73 – -0.31 | 0.005 |
| Zdoy | -0.38 | -0.70 – -0.07 | 0.017 |
| Sex [M] × Year [2010] | 0.99 | 0.13 – 1.86 | 0.024 |
| Sex [M] × Year [2011] | 1.44 | 0.54 – 2.35 | 0.002 |
| Sex [M] × Zdoy | 0.42 | 0.00 – 0.83 | 0.047 |
| Random Effects | |||
| σ2 | 3.04 | ||
| τ00 FishID | 0.38 | ||
| τ11 FishID.Zdoy | 0.16 | ||
| ρ01 FishID | 0.44 | ||
| ICC | 0.15 | ||
| N FishID | 36 | ||
| Observations | 813 | ||
| Marginal R2 / Conditional R2 | 0.095 / 0.231 | ||
simulationOutput <- simulateResiduals(fittedModel=mM, plot=FALSE)
testTemporalAutocorrelation(simulationOutput, time=mMdat$datetime)## Warning: Unknown or uninitialised column: `datetime`.## DHARMa::testTemporalAutocorrelation - no time argument provided, using random times for each data point
##
## Durbin-Watson test
##
## data: simulationOutput$scaledResiduals ~ 1
## DW = 1.9293, p-value = 0.3131
## alternative hypothesis: true autocorrelation is not 0mMreduced <- lmer(log(distday) ~ Sex*Year + Sex*Zdoy + (1|FishID), data=mMdat)
anova(mMreduced, mM)## refitting model(s) with ML (instead of REML)## Data: mMdat
## Models:
## mMreduced: log(distday) ~ Sex * Year + Sex * Zdoy + (1 | FishID)
## mM: log(distday) ~ Sex * Year + Sex * Zdoy + (Zdoy | FishID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## mMreduced 10 3301.2 3348.2 -1640.6 3281.2
## mM 12 3292.8 3349.2 -1634.4 3268.8 12.416 2 0.002013 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Change point models
Change point model outputs are produced in a separate script, and summaries of their results are uploaded here for presentation and for use in the home range analysis.
mcp_summary %>% arrange(Usedmod)## # A tibble: 36 × 18
## i FishID Topmod Usedmod n class Frequency Sex TLmm Wkg
## <dbl> <dbl> <dbl> <dbl> <dbl> <chr> <dbl> <chr> <dbl> <dbl>
## 1 6 20 1 0 5 Stationary 150. F 704 3.1
## 2 9 23 2 0 15 Stationary 150. M 576 1.5
## 3 28 43 3 0 34 Stationary 150. M 616 1.95
## 4 29 44 0 0 31 Stationary 150. M 581 1.75
## 5 31 47 0 0 33 Stationary 150. M 570 1.9
## 6 34 32 0 0 8 Stationary 150. F 705 3.35
## 7 13 28 1 1 33 Migration 150. F 676 2.4
## 8 1 12 2 2 24 Migration 150. M 545 1.38
## 9 5 19 2 2 7 Migration 150. F 674 3
## 10 7 21 2 2 27 Migration 150. F 646 1.95
## # ℹ 26 more rows
## # ℹ 8 more variables: ReleaseDate <dttm>, ReleaseSite <chr>,
## # SurgeryTime.min <chr>, Surgeon <chr>, BatteryEnd <dttm>, Omit <dbl>,
## # Notes <chr>, id <chr>Reorganize fixed effect coefficients from change point models.
# concatenate seasonal change points and filter out non-seasonal models
mcp_fixefs <- lapply(1:length(mcp_topmod), function(x){
summary(mcp_topmod[[x]]$mcmc_post)$quantiles %>% data.frame() %>%
rename(median=X50.) %>%
mutate(name=row.names(.),
id=names(mcp_topmod[x])) %>%
filter(!grepl("\\[", name))
}) %>% bind_rows()
# extract random effects - median change point by year by fish
mcp_cp1yr <- lapply(1:length(mcp_topmod), function(x){
summary(mcp_topmod[[x]]$mcmc_post)$quantiles %>% data.frame() %>%
rename(median=X50.) %>%
mutate(name=row.names(.),
id=names(mcp_topmod[x])) %>%
separate("name", into=c("name", "year"), sep="\\[|\\]", extra="drop") %>%
filter(name=="cp_1"|name=="cp_1_year") %>%
select(id, name, year, median) %>%
mutate(cp1day=median+pull(filter(., name=="cp_1"), median)) %>%
filter(!is.na(year))
}) %>% bind_rows() ## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 6 rows [1, 2, 6, 7, 8,
## 9].## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 9 rows [1, 2, 6, 7, 11,
## 12, 13, 14, 15].
## Expected 2 pieces. Missing pieces filled with `NA` in 9 rows [1, 2, 6, 7, 11,
## 12, 13, 14, 15].
## Expected 2 pieces. Missing pieces filled with `NA` in 9 rows [1, 2, 6, 7, 11,
## 12, 13, 14, 15].## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 6 rows [1, 2, 5, 6, 7,
## 8].
## Expected 2 pieces. Missing pieces filled with `NA` in 6 rows [1, 2, 5, 6, 7,
## 8].## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 6 rows [1, 2, 6, 7, 8,
## 9].## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 9 rows [1, 2, 6, 7, 11,
## 12, 13, 14, 15].## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 6 rows [1, 2, 6, 7, 8,
## 9].## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 9 rows [1, 2, 6, 7, 11,
## 12, 13, 14, 15].## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 6 rows [1, 2, 6, 7, 8,
## 9].## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 9 rows [1, 2, 6, 7, 11,
## 12, 13, 14, 15].## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 6 rows [1, 2, 6, 7, 8,
## 9].
## Expected 2 pieces. Missing pieces filled with `NA` in 6 rows [1, 2, 6, 7, 8,
## 9].## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 9 rows [1, 2, 6, 7, 11,
## 12, 13, 14, 15].
## Expected 2 pieces. Missing pieces filled with `NA` in 9 rows [1, 2, 6, 7, 11,
## 12, 13, 14, 15].
## Expected 2 pieces. Missing pieces filled with `NA` in 9 rows [1, 2, 6, 7, 11,
## 12, 13, 14, 15].## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 6 rows [1, 2, 6, 7, 8,
## 9].## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 9 rows [1, 2, 6, 7, 11,
## 12, 13, 14, 15].## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 6 rows [1, 2, 6, 7, 8,
## 9].## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 9 rows [1, 2, 6, 7, 11,
## 12, 13, 14, 15].## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 6 rows [1, 2, 6, 7, 8,
## 9].## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 9 rows [1, 2, 6, 7, 11,
## 12, 13, 14, 15].
## Expected 2 pieces. Missing pieces filled with `NA` in 9 rows [1, 2, 6, 7, 11,
## 12, 13, 14, 15].
## Expected 2 pieces. Missing pieces filled with `NA` in 9 rows [1, 2, 6, 7, 11,
## 12, 13, 14, 15].
## Expected 2 pieces. Missing pieces filled with `NA` in 9 rows [1, 2, 6, 7, 11,
## 12, 13, 14, 15].## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 6 rows [1, 2, 5, 6, 7,
## 8].## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 9 rows [1, 2, 6, 7, 11,
## 12, 13, 14, 15].## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 2 rows [1, 2].## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 6 rows [1, 2, 6, 7, 8,
## 9].## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 2 rows [1, 2].## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 9 rows [1, 2, 6, 7, 11,
## 12, 13, 14, 15].
## Expected 2 pieces. Missing pieces filled with `NA` in 9 rows [1, 2, 6, 7, 11,
## 12, 13, 14, 15].## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 2 rows [1, 2].## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 6 rows [1, 2, 6, 7, 8,
## 9].## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 9 rows [1, 2, 6, 7, 11,
## 12, 13, 14, 15].mcp_fixefs %>% left_join(mcp_summary) %>%
filter(name=="cp_1"& class=="Migration") %>%
summarize(mean=mean(median), sd=sd(median), min=min(median), max=max(median),
n=n())## Joining with `by = join_by(id)`## mean sd min max n
## 1 162.1114 20.92017 118.5156 193.3826 30How many fish with change point models supported, and do their change
point dates differ by sex? The top models for some fish contained change
point estimates but were rejected as still not supporting migratory
behavior due to lack of fit. Because of this, we need cp_1
parameters from fish tracks classified as “Migration”.
mcp_fixefs %>% left_join(mcp_summary) %>%
filter(name=="cp_1"& class=="Migration") %>%
group_by(Sex) %>%
summarize(mean=mean(median), sd=sd(median), min=min(median), max=max(median),
n=n())## Joining with `by = join_by(id)`## # A tibble: 2 × 6
## Sex mean sd min max n
## <chr> <dbl> <dbl> <dbl> <dbl> <int>
## 1 F 160. 20.8 128. 193. 16
## 2 M 164. 21.6 119. 193. 14mcp_fixefs %>% left_join(mcp_summary) %>%
filter(name=="cp_1"& class=="Migration") %>%
{t.test(.$median ~ .$Sex, var.equal = FALSE)}## Joining with `by = join_by(id)`##
## Welch Two Sample t-test
##
## data: .$median by .$Sex
## t = -0.52384, df = 27.19, p-value = 0.6046
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
## -20.01236 11.86997
## sample estimates:
## mean in group F mean in group M
## 160.2116 164.2827How many of these fish had two change points between (to and from)
two intervals? We want to compare cp_2; just extracting
that parameter is enough to filter to the target subset.
mcp_fixefs %>% left_join(mcp_summary) %>%
filter(name=="cp_2") %>%
summarize(mean=mean(median), sd=sd(median), min=min(median), max=max(median),
n=n())## Joining with `by = join_by(id)`## mean sd min max n
## 1 263.5515 20.07628 211.2955 289.8022 19mcp_fixefs %>% left_join(mcp_summary) %>%
filter(name=="cp_2") %>%
group_by(Sex) %>%
summarize(mean=mean(median), sd=sd(median), min=min(median), max=max(median),
n=n())## Joining with `by = join_by(id)`## # A tibble: 2 × 6
## Sex mean sd min max n
## <chr> <dbl> <dbl> <dbl> <dbl> <int>
## 1 F 254. 22.4 211. 280. 9
## 2 M 272. 13.8 244. 290. 10mcp_fixefs %>% left_join(mcp_summary) %>%
filter(name=="cp_2") %>%
{t.test(.$median ~ .$Sex, var.equal = FALSE)}## Joining with `by = join_by(id)`##
## Welch Two Sample t-test
##
## data: .$median by .$Sex
## t = -2.0673, df = 13.032, p-value = 0.05917
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
## -36.5923290 0.8009524
## sample estimates:
## mean in group F mean in group M
## 254.1327 272.0284How many of the “seasonal” fish tracks had variance-only changes
points. We can look at this by extracting and summarizing the number of
fish for which Usedmod equals 1.
mcp_fixefs %>% left_join(mcp_summary) %>%
filter(name=="cp_1"&Usedmod==1) %>%
group_by(Sex) %>%
summarize(mean=mean(median), sd=sd(median), min=min(median), max=max(median),
n=n())## Joining with `by = join_by(id)`## # A tibble: 1 × 6
## Sex mean sd min max n
## <chr> <dbl> <dbl> <dbl> <dbl> <int>
## 1 F 183. NA 183. 183. 1Here are results for fish with single change point models.
- Overall:
mcp_fixefs %>% left_join(mcp_summary) %>%
filter(name=="cp_1"&Usedmod==2) %>%
summarize(mean=mean(median), sd=sd(median), min=min(median), max=max(median),
n=n())## Joining with `by = join_by(id)`## mean sd min max n
## 1 158.9763 24.76114 118.5156 193.0988 10- By sex:
mcp_fixefs %>% left_join(mcp_summary) %>%
filter(name=="cp_1"&Usedmod==2) %>%
group_by(Sex) %>%
summarize(mean=mean(median), sd=sd(median), min=min(median), max=max(median),
n=n())## Joining with `by = join_by(id)`## # A tibble: 2 × 6
## Sex mean sd min max n
## <chr> <dbl> <dbl> <dbl> <dbl> <int>
## 1 F 163. 23.2 133. 193. 5
## 2 M 155. 28.3 119. 188. 5Seasons
We delineated season using the change point analysis and used these delineations to classify fish locations.
mcp_sf_lines <- lapply(mcp_sf2, function(x) {
summarize(x, do_union=FALSE) %>% st_cast("LINESTRING")
})
pview <- c(407400.0, 4777456.9, 427398.1, 4789913.8 )
p2fun <- function(i){
name <- paste(names(mcp_sf2)[[i]])
sex <- filter(fishids, FishID==as.numeric(names(mcp_sf2)[[i]])) %>% pull(Sex)
tl <- filter(fishids, FishID==as.numeric(names(mcp_sf2)[[i]])) %>% pull(TLmm)
p1 <- mcp_sf2[[i]] %>%
ggplot() +
geom_tile(data=filter(lrw_df, !is.na(value)), aes(x=x, y=y),
fill='gray90') +
geom_sf(data=olb, color="light blue", fill=NA) +
geom_sf(data=mcp_sf_lines[[i]], color='gray70') +
geom_sf(data=slice(mcp_sf2[[i]], which(t_==max(t_)|t_==min(t_))),
shape=c(1, 2), size=5) +
geom_sf(aes(color=doy, group=year)) +
scale_color_viridis_c(name="DOY") +
theme_void() +
coord_sf(xlim=pview[c(1,3)], ylim=pview[c(2,4)])
p2 <- mcp_sf2[[i]] %>%
mutate(season = if (exists('season', where=.)) season else 5) %>%
ggplot(aes(x=t_, y=nsd)) +
ggtitle("") +
geom_vline(xintercept=as.POSIXct(c("2009-01-01 00:00:00",
"2010-01-01 00:00:00",
"2011-01-01 00:00:00",
"2012-01-01 00:00:00")),
lty=2, color="dark gray") +
geom_line(aes(group=factor(year(t_)))) +
geom_point(aes(fill=factor(season)), color="grey50", shape=21, size=3) +
scale_fill_grey(name="Season", na.value="white") +
#scale_color_discrete(name="Year") +
ylim(0, 350) +
coord_cartesian(xlim=as.POSIXct(c("2009-01-01 00:00:00",
"2012-01-01 00:00:00"))) +
ylab(expression(NSD~(km^2))) + xlab("Time") +
theme_classic() +
theme(legend.position="none")
p2marg <- ggMarginal(p2, type="histogram", margin="y")
p <- (p1+p2marg)
return(p)
}
p.seasonalM <- p2fun(12) #FishID #27 fit2, seasonal male
p.seasonalF <- p2fun(23) #FishID #38 fit3, seasonal female
p.seasonalVar <- p2fun(13) #FishID #28 fit1, var-seasonal female
p.aseasonal <- p2fun(29) #FishID #44 Stationary male
fig5 <- p.seasonalM[[1]] + p.seasonalM[[2]] +
p.seasonalF[[1]] + p.seasonalF[[2]] +
p.seasonalVar[[1]] + p.seasonalVar[[2]] +
p.aseasonal[[1]] + p.aseasonal[[2]] +
plot_layout(nrow=4, ncol=2) +
plot_annotation(tag_levels='a', tag_suffix = ')')
ggsave("fig5_pCPM.png", fig5, width=8, height=8, units="in", dpi=600, bg="white")
knitr::include_graphics("fig5_pCPM.png")
pcpdist <- filter(mcp_fixefs, name=="cp_1") %>% left_join(fishids) %>%
ggplot(aes(factor(Sex), median, color=TLmm)) +
geom_violin(fill="white") +
ylab("Change point (day of year)") +
xlab("Sex") +
ggbeeswarm::geom_beeswarm(size=3) +
scale_color_viridis_c() + theme_minimal() + labs(fill=c("P(DiffInt)"))## Joining with `by = join_by(id)`pcpdist## Warning: The following aesthetics were dropped during statistical transformation: colour
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
## the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
## variable into a factor?
Define separate tracks for separate seasons.
# all relocations
dat_hr <- mcp_sf2_tab %>% left_join(fishids) %>% left_join(mcp_summary) %>%
mutate(x_=st_coordinates(.)[,1], y_=st_coordinates(.)[,2]) %>%
st_drop_geometry() %>% as_tibble() %>%
select(x_, y_, t_, sex=Sex, tl=TLmm, kg=Wkg, id, rel=ReleaseSite,
t_rel=ReleaseDate, surgeon=Surgeon, year, season) ## Joining with `by = join_by(id)`
## Joining with `by = join_by(id, FishID, Frequency, Sex, TLmm, Wkg, ReleaseDate,
## ReleaseSite, SurgeryTime.min, Surgeon, BatteryEnd, Omit, Notes)`# all seasons
dat_all <- select(dat_hr, -season) %>%
make_track(x_, y_, t_, id, sex, tl, kg, rel, crs=crs_utm18n) %>%
nest(data =c(x_, y_, t_))
# spring & fall seasons separately
dat_season <- filter(dat_hr, season%in%c(0,1)) %>%
make_track(x_, y_, t_, id, sex, tl, kg, rel, season, crs=crs_utm18n) %>%
nest(data =c(x_, y_, t_))Home Ranges
Organize relocation data into tidy framework employed by package
amt. To do so, condition into a movement track list - a
tibble with one row per fish and a list data column containing
relocations. Next, add a new list column that contains the home-range
estimate for each animal. Estimate home ranges using 4 methods as Signer
et al. (2015) did, truncating the dataset to only tracks with greater
than 9 relocations.
hrlevel <- 0.8
# Home ranges for full relocation dataset
hr_all <- dat_all %>%
mutate(n_relocs=map_dbl(data, ~nrow(.))) %>% filter(n_relocs>9) %>%
mutate(
hr_mcp = map(data, hr_mcp, levels=hrlevel),
hr_locoh = map(data, possibly(~ hr_locoh(., n=ceiling(sqrt(nrow(.))),
levels=hrlevel), NULL)),
hr_kde = map(data, ~ hr_kde(., trast=make_trast(., factor=1.75),
levels=hrlevel)),
hr_akde = map(data, ~ hr_akde(., fit_ctmm(., "auto"), levels=hrlevel))
)
#save(hr_all, file="hr_all.Rdata")
cat(lu(hr_all$id), "fish with all-season home range estimates\n")## 32 fish with all-season home range estimates# Home ranges for tracks split by MCP results
hr_season <- dat_season %>%
mutate(n_relocs=map_dbl(data, ~nrow(.))) %>% filter(n_relocs>9) %>%
mutate(
hr_mcp = map(data, hr_mcp, levels=hrlevel),
hr_locoh = map(data, possibly(~ hr_locoh(., n=ceiling(sqrt(nrow(.))),
levels=hrlevel), NULL)),
hr_kde = map(data, ~ hr_kde(., trast=make_trast(., factor=1.75),
levels=hrlevel)),
hr_akde = map(data, ~ hr_akde(., fit_ctmm(., "auto"), levels=hrlevel))
)
#save(hr_season, file="hr_season.Rdata")
cat(lu(hr_season$id), "fish with season-specific home range estimates\n")## 24 fish with season-specific home range estimatescat(filter(hr_season, season==0) %>% nrow(), "spring home range estimates\n")## 14 spring home range estimatescat(filter(hr_season, season==1) %>% nrow(), "summer home range estimates\n")## 23 summer home range estimatescat(group_by(hr_season, id) %>% summarize(n=n()) %>%
filter(n==2) %>% nrow(), "with both seasons' home range estimates")## 13 with both seasons' home range estimatesHome range models have now been fit and stored in table objects. Pivot to long format for analysis of home range sizes by estimator.
# pivot to long format
hr_all_1 <- hr_all %>% select(-data) %>%
pivot_longer(hr_mcp:hr_akde, names_to="estimator", values_to="hr")
hr_season_1 <- hr_season %>% select(-data) %>%
pivot_longer(hr_mcp:hr_akde, names_to="estimator", values_to="hr")
# Functions to calculate areas (in square m)
do_iso <- function(x) hr_isopleths(x) %>% filter(what=="estimate")
clip_iso <- function(x) st_intersection(x, lrwpoly)
# Calculate home range areas
hr_all_1.area <- hr_all_1 %>%
mutate(hr_area=map(hr, possibly(hr_area, NULL)),
hr_isopleth=map(hr, possibly(do_iso, NULL)),
hr_isopleth_clp=map(hr_isopleth, possibly(clip_iso, NULL))
)
hr_season_1.area <- hr_season_1%>%
mutate(hr_area=map(hr, possibly(hr_area, NULL)),
hr_isopleth=map(hr, possibly(do_iso, NULL)),
hr_isopleth_clp=map(hr_isopleth, possibly(clip_iso, NULL))
)Map home ranges
Use the home ranges from aKDE estimation for visualization.
all_akde <- hr_all_1.area %>% filter(estimator=="hr_akde") %>%
select(-hr_isopleth) %>% unnest(hr_isopleth_clp) %>% st_as_sf()
sp_akde <- hr_season_1.area %>% filter(season==0&estimator=="hr_akde") %>%
select(-hr_isopleth) %>% unnest(hr_isopleth_clp) %>% st_as_sf()
su_akde <- hr_season_1.area %>% filter(season==1&estimator=="hr_akde") %>%
select(-hr_isopleth) %>% unnest(hr_isopleth_clp) %>% st_as_sf()#main plot
plotlocations_f <- ggplot() +
geom_tile(data=rw_df %>% filter(!is.na(value)),
aes(x=x, y=y), fill='gray90', color='gray90') +
geom_tile(data=bathy_df %>% filter(!is.na(value)),
aes(x=x, y=y, fill=value, color=value)) +
scale_fill_viridis_c(name='Depth (m)', direction=-1, end=0.6) +
scale_color_viridis_c(name='Depth (m)', direction=-1, end=0.6) +
theme_void()
# lines for raw HRs
sp_akde_unclip <- hr_season_1.area %>% filter(season==0&estimator=="hr_akde") %>%
select(-hr_isopleth_clp) %>% unnest(hr_isopleth) %>% st_as_sf() %>%
arrange(-area)
su_akde_unclip <- hr_season_1.area %>% filter(season==1&estimator=="hr_akde") %>%
select(-hr_isopleth_clp) %>% unnest(hr_isopleth) %>% st_as_sf() %>%
arrange(-area)
# plot fill for clipped but outlines for raw HRs
p_sp_akde_f <- plotlocations_f +
ggtitle("Spring") +
geom_sf(data=sp_akde_unclip, alpha=0.2, fill="gold", color=NA) +
geom_sf(data=sp_akde_unclip, alpha=0.2, fill=NA, color="white") +
geom_tile(data=lrw_mask_df %>% filter(!is.na(value)),
aes(x=x, y=y), fill="white", color="white") +
geom_sf(data=st_centroid(sp_akde_unclip), shape=3, color="gray20") +
coord_sf(xlim=st_bbox(olb)[c(1,3)], ylim=st_bbox(olb)[c(2,4)]) ## Warning: st_centroid assumes attributes are constant over geometriesp_su_akde_f <- plotlocations_f +
ggtitle("Summer") +
geom_sf(data=su_akde, alpha=0.2, fill="gold", color=NA) +
geom_sf(data=su_akde, alpha=0.2, fill=NA, color="white") +
geom_sf(data=st_centroid(su_akde_unclip), shape=3, color="gray20") +
coord_sf(xlim=st_bbox(olb)[c(1,3)], ylim=st_bbox(olb)[c(2,4)]) +
ggspatial::annotation_north_arrow(which_north="true",
pad_y = unit(0.635, "cm"),
height = unit(1., "cm"),
width = unit(1., "cm"),
style = ggspatial::north_arrow_fancy_orienteering) +
ggspatial::annotation_scale(location = "bl", width_hint = 0.5) ## Warning: st_centroid assumes attributes are constant over geometries#build plot
combined <- p_sp_akde_f / p_su_akde_f + plot_layout(guides="collect")
ggsave("fig6_pHRS.png", combined, width=6.5, height=4.3, dpi=600, bg="white")
knitr::include_graphics("fig6_pHRS.png")
Plot home range sizes
Plot home range size by estimator and sex for comparison. Home ranges are expressed in meters because that is the basic unit of the coordinate reference system, but here we convert to square km.
Aseasonal home range
# function to derive home range sizes
fci <- function(hr.area){
hr.area1 <- hr.area %>%
unnest(cols=hr_area) %>%
mutate(area = area / 1e6)
hr.area1$est_lab<-factor(hr.area1$estimator,
levels=c("hr_mcp", "hr_locoh", "hr_kde", "hr_akde"),
labels=c("MCP", "LoCoH", "KDE", "aKDE"))
#remove list column in order to pivot wider and separate estimates from ci
hr.area2 <- hr.area1 %>% select(-hr) %>%
pivot_wider(names_from=what, values_from=area)
#calculate upper and lower confidence intervals on log scale
ci2 <- hr.area2 %>% group_by(est_lab, sex) %>%
summarise(m = mean(log(estimate)),
sd = sd(log(estimate)),
se = sd(log(estimate)) / sqrt(n()),
me = qt(0.975, n() - 1) * se,
lci = m - me,
uci = m + me,
n=n())
return(list(hr.area=hr.area2, ci.area=ci2))
}
# function to plot home range sizes
fplot.range <- function(hr.area, ci.area){
# gj edits to p1
p1.2 <- hr.area %>%
ggplot(aes(est_lab, estimate, shape = sex)) +
scale_shape_manual(values=c(2, 5)) +
geom_pointrange(aes(x = est_lab, y = estimate, ymin = `lci (0.95)`,
ymax = `uci (0.95)`, shape = sex), inherit.aes = FALSE,
position = position_dodge2(width = 0.5)) +
theme_light() +
theme(axis.title.x = element_blank()) +
labs(y = expression(paste("HRS [", km^2, "]")))
p2.2 <- hr.area %>%
ggplot(aes(est_lab, log(estimate), shape = sex)) +
scale_shape_manual(values=c(2, 5)) +
geom_jitter(alpha = 0.5, position = position_dodge2(width = 0.5)) +
geom_pointrange(aes(x = est_lab, y = m, ymin = lci, ymax = uci, shape=sex),
data = ci.area, inherit.aes = FALSE,
position = position_dodge2(width = 0.5),
size=1) +
theme_light() +
theme(axis.title.x = element_blank()) +
labs(y = expression(paste("HRS [", km^2, "]")))
# need confidence interval of the log response ratio log(meanM/meanF)
# se(log(barm)-log(barf) = sqrt( s^2 /(n*barx)))
hr.logratio <- hr.area %>%
nest(data=-c(est_lab, sex)) %>%
mutate(barx=map_dbl(data, ~mean(log(.$estimate))),
SD=map_dbl(data, ~sd(log(.$estimate))),
n=map_dbl(data, ~nrow(.)),
Sn=(SD^2)/n) %>%
nest(data=-est_lab) %>%
mutate(lr=map_dbl(data, ~ .$barx[.$sex=="M"] - .$barx[.$sex=="F"]),
SElr=map_dbl(data, ~ sqrt(.$Sn[.$sex=="M"] + .$Sn[.$sex=="F"])),
df=map_dbl(data, ~ .$n[.$sex=="M"]+.$n[.$sex=="F"]-2),
ucilr=lr+qt(0.975, df)*SElr,
lcilr=lr-qt(0.975, df)*SElr)
p3.2 <- hr.logratio %>%
ggplot(aes(est_lab, lr)) +
geom_pointrange(aes(ymin=lcilr, ymax=ucilr)) +
geom_hline(yintercept = 0, lty = 2) +
theme_light() +
theme(axis.title.x = element_blank()) +
labs(y = expression(log(HRS[m]) - log(HRS[f])))
(p1.2 / p2.2 / p3.2) + plot_layout(heights = c(0.5, 0.5, 0.5)) +
plot_annotation(tag_levels = "A", tag_suffix = ")")
}
# Run function for non-seasonal home range areas
(hr_all_out<- fci(hr_all_1.area))## `summarise()` has grouped output by 'est_lab'. You can override using the
## `.groups` argument.## $hr.area
## # A tibble: 128 × 14
## id sex tl kg rel n_relocs estimator level hr_isopleth
## <chr> <chr> <dbl> <dbl> <chr> <dbl> <chr> <dbl> <list>
## 1 12 M 545 1.38 Bradburys 24 hr_mcp 0.8 <sf [1 × 4]>
## 2 12 M 545 1.38 Bradburys 24 hr_locoh 0.79 <sf [1 × 4]>
## 3 12 M 545 1.38 Bradburys 24 hr_kde 0.8 <sf [1 × 4]>
## 4 12 M 545 1.38 Bradburys 24 hr_akde 0.8 <sf [1 × 4]>
## 5 14 M 635 1.9 Bradburys 55 hr_mcp 0.8 <sf [1 × 4]>
## 6 14 M 635 1.9 Bradburys 55 hr_locoh 0.85 <sf [1 × 4]>
## 7 14 M 635 1.9 Bradburys 55 hr_kde 0.8 <sf [1 × 4]>
## 8 14 M 635 1.9 Bradburys 55 hr_akde 0.8 <sf [1 × 4]>
## 9 15 F 606 1.8 Bradburys 36 hr_mcp 0.8 <sf [1 × 4]>
## 10 15 F 606 1.8 Bradburys 36 hr_locoh 0.81 <sf [1 × 4]>
## # ℹ 118 more rows
## # ℹ 5 more variables: hr_isopleth_clp <list>, est_lab <fct>, estimate <dbl>,
## # `lci (0.95)` <dbl>, `uci (0.95)` <dbl>
##
## $ci.area
## # A tibble: 8 × 9
## # Groups: est_lab [4]
## est_lab sex m sd se me lci uci n
## <fct> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <int>
## 1 MCP F 1.27 2.64 0.704 1.52 -0.254 2.79 14
## 2 MCP M 0.199 1.98 0.466 0.984 -0.785 1.18 18
## 3 LoCoH F -0.635 1.71 0.457 0.988 -1.62 0.353 14
## 4 LoCoH M -2.53 1.34 0.316 0.667 -3.20 -1.87 18
## 5 KDE F 4.01 2.04 0.545 1.18 2.83 5.19 14
## 6 KDE M 3.11 1.64 0.385 0.813 2.30 3.93 18
## 7 aKDE F 3.74 2.19 0.585 1.26 2.48 5.01 14
## 8 aKDE M 2.70 1.61 0.380 0.801 1.90 3.50 18Spring home range
# Run function for seasonal home range areas split by MCP intervals
(hr_season_out_spring <- hr_season_1.area %>% filter(season==0) %>% fci())## `summarise()` has grouped output by 'est_lab'. You can override using the
## `.groups` argument.## $hr.area
## # A tibble: 56 × 15
## id sex tl kg rel season n_relocs estimator level hr_isopleth
## <chr> <chr> <dbl> <dbl> <chr> <dbl> <dbl> <chr> <dbl> <list>
## 1 12 M 545 1.38 Bradburys 0 11 hr_mcp 0.8 <sf>
## 2 12 M 545 1.38 Bradburys 0 11 hr_locoh 0.82 <sf>
## 3 12 M 545 1.38 Bradburys 0 11 hr_kde 0.8 <sf>
## 4 12 M 545 1.38 Bradburys 0 11 hr_akde 0.8 <sf>
## 5 14 M 635 1.9 Bradburys 0 13 hr_mcp 0.8 <sf>
## 6 14 M 635 1.9 Bradburys 0 13 hr_locoh 0.77 <sf>
## 7 14 M 635 1.9 Bradburys 0 13 hr_kde 0.8 <sf>
## 8 14 M 635 1.9 Bradburys 0 13 hr_akde 0.8 <sf>
## 9 15 F 606 1.8 Bradburys 0 16 hr_mcp 0.8 <sf>
## 10 15 F 606 1.8 Bradburys 0 16 hr_locoh 0.81 <sf>
## # ℹ 46 more rows
## # ℹ 5 more variables: hr_isopleth_clp <list>, est_lab <fct>, estimate <dbl>,
## # `lci (0.95)` <dbl>, `uci (0.95)` <dbl>
##
## $ci.area
## # A tibble: 8 × 9
## # Groups: est_lab [4]
## est_lab sex m sd se me lci uci n
## <fct> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <int>
## 1 MCP F -0.292 2.53 1.03 2.66 -2.95 2.37 6
## 2 MCP M -2.19 2.51 0.886 2.09 -4.28 -0.0927 8
## 3 LoCoH F -2.22 2.72 1.11 2.86 -5.07 0.641 6
## 4 LoCoH M -3.50 2.84 1.00 2.37 -5.87 -1.13 8
## 5 KDE F 2.44 2.28 0.932 2.40 0.0460 4.84 6
## 6 KDE M 0.965 1.33 0.469 1.11 -0.144 2.07 8
## 7 aKDE F 1.84 2.48 1.01 2.60 -0.761 4.44 6
## 8 aKDE M 0.318 1.50 0.531 1.26 -0.937 1.57 8Fall home range
# Run function for seasonal home range areas split by MCP intervals
(hr_season_out_su <- hr_season_1.area %>% filter(season==1) %>% fci())## `summarise()` has grouped output by 'est_lab'. You can override using the
## `.groups` argument.## $hr.area
## # A tibble: 92 × 15
## id sex tl kg rel season n_relocs estimator level hr_isopleth
## <chr> <chr> <dbl> <dbl> <chr> <dbl> <dbl> <chr> <dbl> <list>
## 1 29 M 617 2.1 Bradburys 1 26 hr_mcp 0.8 <sf>
## 2 29 M 617 2.1 Bradburys 1 26 hr_locoh 0.77 <sf>
## 3 29 M 617 2.1 Bradburys 1 26 hr_kde 0.8 <sf>
## 4 29 M 617 2.1 Bradburys 1 26 hr_akde 0.8 <sf>
## 5 12 M 545 1.38 Bradburys 1 13 hr_mcp 0.8 <sf>
## 6 12 M 545 1.38 Bradburys 1 13 hr_locoh 0.92 <sf>
## 7 12 M 545 1.38 Bradburys 1 13 hr_kde 0.8 <sf>
## 8 12 M 545 1.38 Bradburys 1 13 hr_akde 0.8 <sf>
## 9 15 F 606 1.8 Bradburys 1 11 hr_mcp 0.8 <sf>
## 10 15 F 606 1.8 Bradburys 1 11 hr_locoh 0.91 <sf>
## # ℹ 82 more rows
## # ℹ 5 more variables: hr_isopleth_clp <list>, est_lab <fct>, estimate <dbl>,
## # `lci (0.95)` <dbl>, `uci (0.95)` <dbl>
##
## $ci.area
## # A tibble: 8 × 9
## # Groups: est_lab [4]
## est_lab sex m sd se me lci uci n
## <fct> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <int>
## 1 MCP F -0.879 2.81 0.812 1.79 -2.67 0.909 12
## 2 MCP M -1.50 2.03 0.613 1.37 -2.87 -0.139 11
## 3 LoCoH F -1.73 2.63 0.759 1.67 -3.40 -0.0604 12
## 4 LoCoH M -2.86 1.56 0.470 1.05 -3.90 -1.81 11
## 5 KDE F 1.57 3.12 0.901 1.98 -0.416 3.55 12
## 6 KDE M 0.590 2.31 0.696 1.55 -0.959 2.14 11
## 7 aKDE F 1.06 3.00 0.866 1.91 -0.843 2.97 12
## 8 aKDE M 0.182 2.21 0.667 1.49 -1.31 1.67 11HRs for appendix
p1.1 <- hr_season_out_spring$hr.area %>%
ggplot(aes(est_lab, estimate, shape = sex)) +
scale_shape_manual(values=c(2, 5)) +
geom_pointrange(aes(x = est_lab, y = estimate, ymin = `lci (0.95)`,
ymax = `uci (0.95)`, shape = sex), inherit.aes = FALSE,
position = position_dodge2(width = 0.5)) +
theme_light() +
theme(axis.title.x = element_blank()) +
labs(y = expression(paste("HRS [", km^2, "]")))
p1.2 <- hr_season_out_su$hr.area %>%
ggplot(aes(est_lab, estimate, shape = sex)) +
scale_shape_manual(values=c(2, 5)) +
geom_pointrange(aes(x = est_lab, y = estimate, ymin = `lci (0.95)`,
ymax = `uci (0.95)`, shape = sex), inherit.aes = FALSE,
position = position_dodge2(width = 0.5)) +
theme_light() +
theme(axis.title.x = element_blank()) +
labs(y = expression(paste("HRS [", km^2, "]")))
p2.1 <- hr_season_out_spring$hr.area %>%
ggplot(aes(est_lab, log(estimate), shape = sex)) +
scale_shape_manual(values=c(2, 5)) +
geom_jitter(alpha = 0.5, position = position_dodge2(width = 0.5)) +
geom_pointrange(aes(x = est_lab, y = m, ymin = lci, ymax = uci, shape=sex),
data = hr_season_out_spring$ci.area, inherit.aes = FALSE,
position = position_dodge2(width = 0.5),
size=1) +
theme_light() +
theme(axis.title.x = element_blank()) +
labs(y = expression(paste("HRS [", km^2, "]")))
p2.2 <- hr_season_out_su$hr.area %>%
ggplot(aes(est_lab, log(estimate), shape = sex)) +
scale_shape_manual(values=c(2, 5)) +
geom_jitter(alpha = 0.5, position = position_dodge2(width = 0.5)) +
geom_pointrange(aes(x = est_lab, y = m, ymin = lci, ymax = uci, shape=sex),
data = hr_season_out_su$ci.area, inherit.aes = FALSE,
position = position_dodge2(width = 0.5),
size=1) +
theme_light() +
theme(axis.title.x = element_blank()) +
labs(y = expression(paste("HRS [", km^2, "]")))
# need confidence interval of the log response ratio log(meanM/meanF)
# se(log(barm)-log(barf) = sqrt( s^2 /(n*barx)))
hr.logratio1 <- hr_season_out_spring$hr.area %>%
nest(data=-c(est_lab, sex)) %>%
mutate(barx=map_dbl(data, ~mean(log(.$estimate))),
SD=map_dbl(data, ~sd(log(.$estimate))),
n=map_dbl(data, ~nrow(.)),
Sn=(SD^2)/n) %>%
nest(data=-est_lab) %>%
mutate(lr=map_dbl(data, ~ .$barx[.$sex=="M"] - .$barx[.$sex=="F"]),
SElr=map_dbl(data, ~ sqrt(.$Sn[.$sex=="M"] + .$Sn[.$sex=="F"])),
df=map_dbl(data, ~.$n[.$sex=="M"]+.$n[.$sex=="F"]-2),
ucilr=lr+qt(0.975, df)*SElr,
lcilr=lr-qt(0.975, df)*SElr)
hr.logratio2 <- hr_season_out_su$hr.area %>%
nest(data=-c(est_lab, sex)) %>%
mutate(barx=map_dbl(data, ~mean(log(.$estimate))),
SD=map_dbl(data, ~sd(log(.$estimate))),
n=map_dbl(data, ~nrow(.)),
Sn=(SD^2)/n) %>%
nest(data=-est_lab) %>%
mutate(lr=map_dbl(data, ~ .$barx[.$sex=="M"] - .$barx[.$sex=="F"]),
SElr=map_dbl(data, ~ sqrt(.$Sn[.$sex=="M"] + .$Sn[.$sex=="F"])),
df=map_dbl(data, ~.$n[.$sex=="M"]+.$n[.$sex=="F"]-2),
ucilr=lr+qt(0.975, df)*SElr,
lcilr=lr-qt(0.975, df)*SElr)
p3.1 <- hr.logratio1 %>%
ggplot(aes(est_lab, lr)) +
geom_pointrange(aes(ymin=lcilr, ymax=ucilr)) +
geom_hline(yintercept = 0, lty = 2) +
theme_light() +
theme(axis.title.x = element_blank()) +
labs(y = expression(log(HRS[m]) - log(HRS[f])))
p3.2 <- hr.logratio2 %>%
ggplot(aes(est_lab, lr)) +
geom_pointrange(aes(ymin=lcilr, ymax=ucilr)) +
geom_hline(yintercept = 0, lty = 2) +
theme_light() +
theme(axis.title.x = element_blank()) +
labs(y = expression(log(HRS[m]) - log(HRS[f])))
# build plot
pspring <- (p1.1 + p2.1 + p3.1) +
plot_layout(ncol=1)
psummer <- (p1.2 + p2.2 + p3.2) +
plot_layout(ncol=1)
pcombo <- (pspring | psummer) +
plot_layout(ncol=2) +
plot_annotation(tag_levels = "A", tag_suffix = ")")
ggsave("figS4_HRs.png", pcombo, width=9, height=5.3, units="in", dpi=600,
bg="white")## Warning: Removed 42 rows containing missing values (`geom_segment()`).## Warning: Removed 69 rows containing missing values (`geom_segment()`).knitr::include_graphics("figS4_HRs.png")
Aggregate and summarize HR estimates
Create a table summarizing home range size esitmates by season, sex, and estimator
(hr.summary.table <- bind_rows(mutate(hr_all_out$hr.area, season=99),
hr_season_out_spring$hr.area,
hr_season_out_su$hr.area) %>%
mutate(sex=factor(sex),
season=factor(season),
estimator=factor(estimator)) %>%
dplyr::select(season, id:level, est_lab, estimate) %>%
group_by(sex, estimator, season) %>%
summarise(m = mean(log(estimate)),
sd = sd(log(estimate)),
se = sd(log(estimate)) / sqrt(n()),
me = qt(0.975, n() - 1) * se,
lci = m - me,
uci = m + me,
n=n()) %>%
mutate(hat_hr=exp(m), hat_lci=exp(lci), hat_uci=exp(uci)) %>%
arrange(estimator, season, sex))## `summarise()` has grouped output by 'sex', 'estimator'. You can override using
## the `.groups` argument.## # A tibble: 24 × 13
## # Groups: sex, estimator [8]
## sex estimator season m sd se me lci uci n hat_hr
## <fct> <fct> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <int> <dbl>
## 1 F hr_akde 0 1.84 2.48 1.01 2.60 -0.761 4.44 6 6.28
## 2 M hr_akde 0 0.318 1.50 0.531 1.26 -0.937 1.57 8 1.37
## 3 F hr_akde 1 1.06 3.00 0.866 1.91 -0.843 2.97 12 2.90
## 4 M hr_akde 1 0.182 2.21 0.667 1.49 -1.31 1.67 11 1.20
## 5 F hr_akde 99 3.74 2.19 0.585 1.26 2.48 5.01 14 42.1
## 6 M hr_akde 99 2.70 1.61 0.380 0.801 1.90 3.50 18 14.9
## 7 F hr_kde 0 2.44 2.28 0.932 2.40 0.0460 4.84 6 11.5
## 8 M hr_kde 0 0.965 1.33 0.469 1.11 -0.144 2.07 8 2.63
## 9 F hr_kde 1 1.57 3.12 0.901 1.98 -0.416 3.55 12 4.79
## 10 M hr_kde 1 0.590 2.31 0.696 1.55 -0.959 2.14 11 1.80
## # ℹ 14 more rows
## # ℹ 2 more variables: hat_lci <dbl>, hat_uci <dbl>hr.summary.table %>%
mutate(cell=paste0(round(hat_hr, 3), "(",
round(hat_lci, 3), " - ",
round(hat_uci, 3), ")")) %>%
dplyr::select(sex, estimator, season, cell) %>%
pivot_wider(names_from=c(season, sex), values_from=cell) %>%
kableExtra::kable()| estimator | 0_F | 0_M | 1_F | 1_M | 99_F | 99_M |
|---|---|---|---|---|---|---|
| hr_akde | 6.282(0.467 - 84.445) | 1.375(0.392 - 4.826) | 2.897(0.43 - 19.495) | 1.199(0.271 - 5.306) | 42.121(11.891 - 149.195) | 14.938(6.706 - 33.277) |
| hr_kde | 11.499(1.047 - 126.279) | 2.626(0.866 - 7.961) | 4.787(0.66 - 34.739) | 1.805(0.383 - 8.502) | 55.238(17.023 - 179.24) | 22.509(9.982 - 50.756) |
| hr_locoh | 0.109(0.006 - 1.899) | 0.03(0.003 - 0.324) | 0.177(0.033 - 0.941) | 0.057(0.02 - 0.164) | 0.53(0.197 - 1.423) | 0.079(0.041 - 0.155) |
| hr_mcp | 0.747(0.052 - 10.672) | 0.112(0.014 - 0.911) | 0.415(0.069 - 2.481) | 0.222(0.057 - 0.871) | 3.552(0.776 - 16.266) | 1.22(0.456 - 3.265) |
Home range size analysis
hr2.area.season <- bind_rows(hr_season_out_spring$hr.area,
hr_season_out_su$hr.area) %>%
mutate(sex=factor(sex),
season=factor(season),
estimator=factor(estimator)) %>%
dplyr::select(id:level, est_lab, estimate)
m.hr <- lmer(log(estimate) ~ sex*season + (1|est_lab), data=hr2.area.season)anova(m.hr, ddf="Satterthwaite")## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## sex 51.361 51.361 1 141 9.0508 0.003111 **
## season 0.471 0.471 1 141 0.0830 0.773637
## sex:season 3.537 3.537 1 141 0.6233 0.431163
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1margdf <- emmip(m.hr, ~ sex|season, lmer.df="satterthwaite", plotit=FALSE)
pd <- position_dodge(.3)
p <- ggplot(data=margdf, aes(x=season, y=exp(yvar), fill=sex, group=sex)) +
geom_errorbar(aes(ymin=exp(yvar-SE), ymax=exp(yvar+SE)),
color='black', width=.2, position=pd) +
geom_line(position=pd) +
geom_point(color='black', size=8, position=pd, shape=21) +
scale_fill_grey() +
ylab(expression("Home Range Size "~(km^2))) +
scale_y_log10() +
scale_x_discrete(labels=c("0"="Spring", "1"="Summer")) +
theme_classic() +
guides(shape=guide_legend(title="Sex")) +
theme(axis.title=element_text(size = rel(1.5)),
axis.text=element_text(size = rel(1.5), color="black"),
legend.title=element_text(size = rel(1.3)),
legend.text=element_text(size = rel(1.3)),
axis.title.x=element_blank())
ggsave("fig7_Ameans.png", p, width=6.5, height=4, units="in", dpi=600,
bg="white")
knitr::include_graphics("fig7_Ameans.png")
A tabular summary of regression model output for the appendix.
tab_model(m.hr, df.method="satterthwaite")| log(estimate) | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | 0.44 | -2.16 – 3.05 | 0.675 |
| sex [M] | -1.54 | -2.82 – -0.27 | 0.018 |
| season [1] | -0.44 | -1.62 – 0.74 | 0.463 |
| sex [M] × season [1] | 0.64 | -0.97 – 2.25 | 0.431 |
| Random Effects | |||
| σ2 | 5.67 | ||
| τ00 est_lab | 3.00 | ||
| ICC | 0.35 | ||
| N est_lab | 4 | ||
| Observations | 148 | ||
| Marginal R2 / Conditional R2 | 0.039 / 0.371 | ||
Summer resorts
First, establish a grid over which we may produce count summaries of relocations per fish per season per year. We specified a hexagonal-cell grid (area ~2.5 km per cell) over Oneida Lake and retained all cells in which any fish was observed during our study. Results of seasonal home range analysis suggested that a 2.5 km^2 grid cell could contain the average home range size in spring or summer for males or females.
# create hexagonal grid over Oneida Lake detection locations
g <- hr_season %>% unnest(data) %>%
mutate(year=year(t_), yearf=factor(year)) %>%
dplyr::select(id, sex, season, year, yearf, x_:t_) %>%
st_as_sf(coords=c(y="x_", x="y_"), remove=FALSE, crs=crs_utm18n) %>%
st_make_grid(square=FALSE, n = c(10,10)) %>% st_sf() %>%
filter(st_intersects(., summarize(fishobs) %>% st_cast("MULTIPOINT"),
sparse=FALSE) %>% as.vector()) %>%
mutate(cell=1:nrow(.), area=st_area(.))Starting from the data table hr_season, create a new
table hr_season_summerresort that adds columns for data
list cells subset by year (data09, data10,
data11) and hr_mcp list cells by year
(hr_mcp09, hr_mcp10,
hr_mcp11).
hr_summerresort <- hr_season %>%
mutate(data09=map(data, ~subset(., year(.$t_)==2009)),
n09=map_dbl(data09, ~nrow(.)),
hr_mcp09=ifelse(n09>2, map_dbl(data09, ~nrow(.), levels=1), NA),
data10=map(data, ~subset(., year(.$t_)==2010)),
n10=map_dbl(data10, ~nrow(.)),
hr_mcp10=ifelse(n10>2, map_dbl(data10, ~nrow(.), levels=1), NA),
data11=map(data, ~subset(., year(.$t_)==2011)),
n11=map_dbl(data11, ~nrow(.)),
hr_mcp11=ifelse(n11>2, map_dbl(data11, ~nrow(.), levels=1), NA),
datacell=map(data,
~st_intersection(
st_as_sf(., coords=c(y="x_", x="y_"),
remove=FALSE,
crs=crs_utm18n),
g)
),
cellct=map(datacell,
~st_drop_geometry(.) %>%
mutate(year=year(t_),
yearf=factor(year,
levels=c("2009", "2010", "2011"))) %>%
group_by(year, yearf, cell) %>%
summarize(n=n()) %>%
ungroup()
),
cellp=map(cellct,
~select(., -yearf) %>%
full_join(expand.grid(year=c(2009:2011),
cell=unique(.$cell)),
by=join_by(year, cell)) %>%
replace(is.na(.), 0) %>%
group_by(year) %>%
mutate(p=n/sum(n)) %>% ungroup() %>%
pivot_wider(id_cols=cell,
names_from=year,
values_from=p,
names_prefix='p') %>%
mutate(srpq0910=sqrt(p2009*p2010),
srpq1011=sqrt(p2011*p2010))
),
bc0910=map_dbl(cellp, ~summarize(., sum(srpq0910)) %>% pull()),
bc1011=map_dbl(cellp, ~summarize(., sum(srpq1011)) %>% pull())
) %>%
rowwise() %>%
mutate(bcmean=mean(c(bc0910, bc1011), na.rm=TRUE))## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.## Warning: There were 37 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `datacell = map(...)`.
## Caused by warning:
## ! attribute variables are assumed to be spatially constant throughout all geometries
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 36 remaining warnings.Plot summer locations and cell distributions for individual fish. Ladies first.
scale_color_years <- function(...){
ggplot2:::manual_scale(
'color',
values = setNames(c("#f98e09", "#57106e", "#000004"),
#c("2009", "2010", "2011")),
c("2011", "2010", "2009")),
...
)
}
hr_summerresortM <- hr_summerresort %>%
filter(sex=="M", season==1&!is.na(bcmean)) %>% arrange(-bcmean)
psrlistM <- lapply(1:nrow(hr_summerresortM), function(i){
sra <- hr_summerresortM[i,] %>%
unnest(data) %>%
st_as_sf(coords=c("x_", "y_"), crs=crs_utm18n)
this <- ggplot(sra) +
geom_text(aes(x=411000, y=4790200, label=
paste0("ID = ", hr_summerresortM[i,"id"],
", BC = ", round(hr_summerresortM[i,"bcmean"], 2))),
size=2.5, hjust=0) +
geom_sf(data=g, fill=NA, color="grey90") +
geom_text(aes(x=408000, y=4782000),
label=paste0("2009, n=", sum(year(sra$t_)==2009)),
size=2, hjust=0, parse=FALSE, color="#000004") +
geom_text(aes(x=408000, y=4781000),
label=paste0("2010, n=", sum(year(sra$t_)==2010)),
size=2, hjust=0, parse=FALSE, color="#57106e") +
geom_text(aes(x=408000, y=4780000),
label=paste0("2011, n=", sum(year(sra$t_)==2011)),
size=2, hjust=0, parse=FALSE, color="#f98e09") +
geom_sf(data=olb, fill=NA, color="grey50") +
geom_sf(aes(color=factor(year(t_))), shape=3, size=3, alpha=1) +
scale_color_years(name="Year") +
coord_sf(xlim=st_bbox(g)[c(1,3)], ylim=st_bbox(g)[c(2,4)]) +
theme_void() +
theme(plot.margin=unit(c(0,0,0,0), "cm"),
plot.title=element_text(size=10)) +
theme(legend.position="none")
return(this)
})
for (x in 1:(length(psrlistM)-1)){
psrlistM[[x]] <- psrlistM[[x]] +
guides(shape="none", color="none", fill="none") +
theme(legend.position="none")
}
fig9 <- wrap_plots(psrlistM) + plot_layout(ncol=3, guides="collect") +
plot_annotation(tag_levels="a", tag_suffix=")")
ggsave("fig9_pBCM.png", fig9, width=6.5, height=6, dpi=300, bg="white")
knitr::include_graphics("fig9_pBCM.png")
hr_summerresortF <- hr_summerresort %>%
filter(sex=="F", season==1&!is.na(bcmean)) %>% arrange(-bcmean)
psrlistF <- lapply(1:nrow(hr_summerresortF), function(i){
sra <- hr_summerresortF[i,] %>%
unnest(data) %>%
st_as_sf(coords=c("x_", "y_"), crs=crs_utm18n)
this <- ggplot(sra) +
geom_text(aes(x=411000, y=4790200, label=
paste0("ID = ", hr_summerresortF[i,"id"],
", BC = ", round(hr_summerresortF[i,"bcmean"], 2))),
size=2.5, hjust=0) +
geom_sf(data=g, fill=NA, color="grey90") +
geom_text(aes(x=408000, y=4782000),
label=paste0("2009, n=", sum(year(sra$t_)==2009)),
size=2, hjust=0, parse=FALSE, color="#000004") +
geom_text(aes(x=408000, y=4781000),
label=paste0("2010, n=", sum(year(sra$t_)==2010)),
size=2, hjust=0, parse=FALSE, color="#57106e") +
geom_text(aes(x=408000, y=4780000),
label=paste0("2011, n=", sum(year(sra$t_)==2011)),
size=2, hjust=0, parse=FALSE, color="#f98e09") +
geom_sf(data=olb, fill=NA, color="grey50") +
geom_sf(aes(color=factor(year(t_))), shape=3, size=3, alpha=1) +
scale_color_years(name="Year") +
coord_sf(xlim=st_bbox(g)[c(1,3)], ylim=st_bbox(g)[c(2,4)]) +
theme_void() +
theme(plot.margin=unit(c(0,0,0,0), "cm"),
plot.title=element_text(size=10)) +
theme(legend.position="none")
return(this)
})
for (x in 1:(length(psrlistF)-1)){
psrlistF[[x]] <- psrlistF[[x]] +
guides(shape="none", color="none", fill="none") +
theme(legend.position="none")
}
fig10 <- wrap_plots(psrlistF) + plot_layout(ncol=3, guides="collect") +
plot_annotation(tag_levels="a", tag_suffix=")")
ggsave("fig10_BCF.png", fig10, width=6.5, height=6, dpi=300, bg="white")
knitr::include_graphics("fig10_BCF.png")
fig8 <- hr_summerresort %>% filter(season==1&!is.na(bcmean)) %>%
ggplot(aes(x=sex, y=bcmean)) +
geom_point(aes(fill=sex), size=3, shape=21, color="grey20",
position=position_jitter(width=0.2)) +
geom_boxplot(outlier.color=NA, fill=NA, color="grey20") +
scale_fill_grey(name="Sex") +
ylab("Bhattacharyya's coeffcient (BC)") + xlab("") +
ylim(0,1) +
theme_classic()
ggsave("fig8_pba_overlap.png", fig8, width=3, height=3, dpi=300, bg="white")## Warning: Removed 1 rows containing missing values (`geom_point()`).knitr::include_graphics("fig8_pba_overlap.png")
hr_summerresort %>% filter(season==1&!is.na(bcmean)) %>%
group_by(sex) %>%
summarize(n=n(), mean=mean(bcmean), sd=sd(bcmean),
q25=quantile(bcmean, probs=c(0.25)),
q75=quantile(bcmean, probs=c(0.75)),
median=median(bcmean), min=min(bcmean), max=max(bcmean), )## # A tibble: 2 × 9
## sex n mean sd q25 q75 median min max
## <chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 F 12 0.672 0.294 0.498 0.895 0.752 0.167 1
## 2 M 11 0.693 0.360 0.666 0.939 0.822 0 0.964Sys.time()-t0## Time difference of 6.737476 minstestjags()## You are using R version 4.3.1 (2023-06-16 ucrt) on a windows machine,
## with the RTerm GUI
## JAGS version 4.3.1 found successfully using the command '/Program
## Files/JAGS/JAGS-4.3.1/x64/bin/jags-terminal.exe'
## The rjags package is installedsessionInfo()## R version 4.3.1 (2023-06-16 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19045)
##
## Matrix products: default
##
##
## locale:
## [1] LC_COLLATE=English_United States.utf8
## [2] LC_CTYPE=English_United States.utf8
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.utf8
##
## time zone: America/New_York
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] DHARMa_0.4.6 sjlabelled_1.2.0 sjmisc_2.8.9
## [4] sjPlot_2.8.15 emmeans_1.9.0 lmerTest_3.1-3
## [7] lme4_1.1-35.1 Matrix_1.6-0 patchwork_1.2.0
## [10] runjags_2.2.2-1.1 amt_0.2.1.0 raster_3.6-26
## [13] sp_2.1-2 sf_1.0-14 readxl_1.4.3
## [16] kableExtra_1.3.4.9000 ggpubr_0.6.0 ggExtra_0.10.1
## [19] lubridate_1.9.3 forcats_1.0.0 stringr_1.5.1
## [22] dplyr_1.1.2 purrr_1.0.2 readr_2.1.4
## [25] tidyr_1.3.0 tibble_3.2.1 ggplot2_3.4.4
## [28] tidyverse_2.0.0
##
## loaded via a namespace (and not attached):
## [1] rstudioapi_0.15.0 jsonlite_1.8.8 datawizard_0.9.1
## [4] magrittr_2.0.3 ggbeeswarm_0.7.2 TH.data_1.1-2
## [7] estimability_1.4.1 farver_2.1.1 nloptr_2.0.3
## [10] rmarkdown_2.25 ragg_1.2.7 vctrs_0.6.3
## [13] minqa_1.2.6 effectsize_0.8.6 terra_1.7-65
## [16] rstatix_0.7.2 webshot_0.5.5 htmltools_0.5.7
## [19] broom_1.0.5 cellranger_1.1.0 sass_0.4.8
## [22] KernSmooth_2.23-22 bslib_0.6.1 sandwich_3.1-0
## [25] zoo_1.8-12 cachem_1.0.8 mime_0.12
## [28] lifecycle_1.0.4 pkgconfig_2.0.3 R6_2.5.1
## [31] fastmap_1.1.1 rbibutils_2.2.16 shiny_1.8.0
## [34] digest_0.6.34 numDeriv_2016.8-1.1 colorspace_2.1-0
## [37] RSpectra_0.16-1 textshaping_0.3.7 labeling_0.4.3
## [40] fansi_1.0.4 timechange_0.2.0 httr_1.4.7
## [43] abind_1.4-5 compiler_4.3.1 proxy_0.4-27
## [46] bit64_4.0.5 withr_3.0.0 backports_1.4.1
## [49] carData_3.0-5 DBI_1.2.1 performance_0.10.8
## [52] highr_0.10 ggsignif_0.6.4 MASS_7.3-60.0.1
## [55] sjstats_0.18.2 classInt_0.4-10 loo_2.6.0
## [58] tools_4.3.1 units_0.8-4 vipor_0.4.7
## [61] lmtest_0.9-40 beeswarm_0.4.0 httpuv_1.6.13
## [64] glue_1.6.2 nlme_3.1-162 promises_1.2.1
## [67] grid_4.3.1 checkmate_2.3.1 generics_0.1.3
## [70] gtable_0.3.4 tzdb_0.4.0 class_7.3-22
## [73] hms_1.1.3 xml2_1.3.6 car_3.1-2
## [76] utf8_1.2.3 pillar_1.9.0 ctmm_1.2.0
## [79] vroom_1.6.3 ggspatial_1.1.9 later_1.3.2
## [82] robustbase_0.99-1 splines_4.3.1 lattice_0.22-5
## [85] FNN_1.1.4 bit_4.0.5 survival_3.5-7
## [88] tidyselect_1.2.0 Gmedian_1.2.7 sfheaders_0.4.3
## [91] miniUI_0.1.1.1 knitr_1.45 svglite_2.1.3
## [94] xfun_0.39 expm_0.999-9 matrixStats_1.2.0
## [97] DEoptimR_1.1-3 stringi_1.7.12 yaml_2.3.8
## [100] boot_1.3-28.1 evaluate_0.23 codetools_0.2-19
## [103] cli_3.6.1 parameters_0.21.3 xtable_1.8-4
## [106] systemfonts_1.0.5 Rdpack_2.6 munsell_0.5.0
## [109] jquerylib_0.1.4 modelr_0.1.11 Rcpp_1.0.11
## [112] ggeffects_1.3.4 png_0.1-8 coda_0.19-4
## [115] parallel_4.3.1 ellipsis_0.3.2 mcp_0.3.3
## [118] bayestestR_0.13.1 viridisLite_0.4.2 mvtnorm_1.2-4
## [121] scales_1.3.0 e1071_1.7-13 crayon_1.5.2
## [124] insight_0.19.7 rlang_1.1.1 rvest_1.0.3
## [127] multcomp_1.4-25